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^^^ ■ The global-scale interior magnetic field Bi needed to account for the Sun's observed 

\ difTerential rotation can be effective only if confined below the convection zone in all 

I— !■ latitudes including, most critically, the polar caps. Axisymmetric solutions are obtained 

P^ I to the nonlinear magnetohydrodynamic equations showing that such polar confinement 

^0 ■ can be brought about by a very weak downwelling fiow U ^ lO^^cms"^ over each pole. 

(— I I Such downwelling is consistent with the helioseismic evidence. All three components of 

^ I' the magnetic field B decay exponentially with altitude across a thin, laminar "magnetic 

JL \ confinement layer" located at the bottom of the tachocline and permeated by spiralling 

;_( i field lines. With realistic parameter values, the thickness of the confinement layer ^ 

jj^ ' 10~ of the Sun's radius, the thickness scale being the magnetic advection-diffusion 

C^ i scale S = rj/U where the magnetic (ohmic) diffusivity 77 « 4.1 x lO^cm^s^^. Alongside 

baroclinic effects and stable thermal stratification, the solutions take into account the 
^SJ . stable compositional stratification of the helium settling layer, if present as in today's 

K^ ' Sun, and the small diffusivity of helium through hydrogen, x ~ 0-9 ^ lO^cm^s^^. The 

^^ . small value of x relative to 77 produces a double boundary-layer structure in which a 

^1 ' "helium sublayer" of smaller vertical scale (x/vY^^^ is sandwiched between the top of 

ly-v . the helium settling layer and the rest of the confinement layer. Solutions are obtained 

using both semi-analytical and purely numerical, finite-difference techniques. The con- 
finement-layer fiows are magnetostrophic to excellent approximation. More precisely, the 
^^ ' principal force balances are between Lorentz, Coriolis, pressure-gradient and buoyancy 

forces, with relative accelerations negligible to excellent approximation. Viscous forces 
are also negligible, even in the helium sublayer where shears are greatest. This is despite 
the kinematic viscosity being somewhat greater than x- We discuss how the confinement 
/\ ' layers at each pole might fit into a global dynamical picture of the solar tachocline. That 

picture, in turn, suggests a new insight into the early Sun and into the longstanding 
enigma of solar lithium depletion. 



1. Introduction 

This paper analyses a new family of laminar magnetostrophic flows that may be im- 
portant for confining the interior magnetic field Bi needed to explain the Sun's differen- 
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tial rotation. As illustrated in figure 1(a) the differential rotation observed within the 
convection zone goes over into near-solid rotation within the radiative, stably strati- 
fied interior, via a thin shear layer called the "tachocline" much of whi ch is also stably 



stratified. The need for the interior field Bj has been argued elsewhere (|McIntvrelll994 : 



Riidiger fc Kitchatinovlll997l: iGough fc Mclntvrelll998i hereafter GM98); the main argu- 



ments are briefly recalled below. The observational evidence together with many ideas 
about the tachocline are reviewed and further referenced in a recent major compendium. 
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The Solar Tachocline ( Hughes et al. 2007), and further discussed m the second edition 
of Mestel's Stellar Magnetism (Mestel 2011) . 

By a "confined" Bi we mean a field most if not all of whose lines are contained beneath 
the convection zone, and held there against magnetic (ohmic) diffusion. Such confine- 
ment is well kno wn to be necessary in order for the field to help enforce solid rotation in 



the interior 



Mc 



Ferraro 



19371: iMestel fc Weisslll987l: ICharbonneau fc MacGregoij[T993 



MacGregor &: Charbonneaulll999f ). and thereby keep the tachocline thin. Confinement 



against magnetic diffusion requires fluid motion. So, besides magnetic effects, a realistic 
theory of confinement must take account of Coriolis effects, stable stratification, baro- 
clinicity, and thermal relaxation. Without these effects we cannot correctly describe, for 
instance, the overall torque balance, which necessarily involves mean meridional circula- 
tions (MMCs) as well as Maxwell stresses. 

The first attempt at a tachocli ne theory was that of 'Spiegel fc ZahnI ( 1992 ). It included 
all the above effects except Bi. iRiidiger fc Kitchatinov. (19971 ) included Bi but omitted 
the other effects. The first attemp t to include all of them wa s that of GM98, in a line of 
investigation further developed bv lGaraud fc GaraudI ( 20081 ). Meanwhile, the dynamical 
importance of compositional as well as therma l stratification (e.g. lMestellll953l ) was sug- 
gested for tachocline theories ( Mclntvrd 12007 ). In particular, the helium settling layer 
beneath the tachocline is nearly impermeable to MMCs because of the small diffusivity of 
helium through hydrogen. This n ear-impermeability of compositionally stratified regions 
has been called the "mu-choke" ( Mestel fc Mosall986 ). The reality of the Sun's helium 
settling layer is stron gly indicated both b y standard solar-evolution models and by helio 



seisinologv (e.g. Christensen-Dalsgaard et al. 1993t ICiacio et aLl 119971 : lElliott fc Cough 



1999t IChristensen-Dalsgaard fc Thompson! 120071) Ifl As will be seen, this combination of 



circumstances gives rise to some new and interesting fluid dynamics. 

The need for the interior fleld Bi arises from a well known difficulty with non-magnetic 
theories. They tend to spread the strong differential rotation of the convection zone 
down into the radiative interior. Although sometimes disputed, this is a robust and 
well- understood consequence of thermal relaxation, interacting with Corio l is effects and 
gyroscopically-pumped MMCs (iHavnes et al\ Il99ll: ISpiegel fc ZahnI Il992l ; lElliottllT997 : 
iMcIn tvre 200T; Garaud fc Bru mmellll200 8') . As shown bv lSpiegel fc ZahnI and confirmed 
by lElliott . this downward spreading or burrowing would have produced a tachocline far 
thicker than observed. The accompanying MMC, acting throughout the Sun's lifetime, 
would also have prevented the helium settling layer from forming. 

To counter the burrowing tendency and to allow the interior to rotate solidly, angu- 
lar momentum has to be transported somehow from the low-latitude tachocline to the 
high-latitud e tachocline. Th e non-magnetic horizontal eddy viscosity proposed for this 
purpose by ISpiegel fc ZahnI is inconsistent with the properties of non- magnetic strati- 
fied t urbulence known from many studies of the terrestrial atmosphere (|McIntvrelll994 . 



20031 fc refs.). Angular- momentum transport by internal g ravity waves is a physicall 



possible alternative (e.g. | Schatzmanlll993 : Zahn et al. 1997t iRogers fc Clatzmaier 200 






Charbonnel fc Talonll2007i fc refs.). However, it is highly improbable as the main mech- 
anism because, by itself, it has no natural tende ncy to produce solid rotation at all 
latitudes and depths (e.g. lPlumb fc McEwanlll978l ). 

A suitably-shaped magnetic field can, by contrast, naturally produce the required an- 
gular momentum transport, via the Alfvenic elasticity of the field lines. A suitable shape 
is one in which the field lines link low latitudes to high latitudes within the tachocline. 
The simplest such shape — simplest by virtue of its axisymmetry — is that suggested 



f Also Christensen-Dalsgaard & Gough 2011, in preparation. 
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Figure 1 . (a) The Sun's differential r otation deduced from lielioseismic data using inverse metfi- 
ods (adapted from lScIiou et 00119981 ). Tfie radiative interior rotates approximately solidly with 
angular velocity f2i = 2.7 x 10~^s~^, or 435 nHz. Within the convection zone, the angular ve- 
locity increases with colatitude through 350, 400, 450 nHz (heavy contours) to a maximum just 
under 470 nHz at the equator. 

(b) Schematic illustration showing the top of the radiative interior (inner sphere) and the 
time-averaged magnetic field threading the tachocline just above. The cutaway outer sphere 
indicates the top of the tachocline, whose depth has been exaggerated. Poloidal magnetic field 
lines emerge from the interior in high latitudes and are wound up into their curved shapes by the 
tachocline's differential rotation, acting against turbulent eddy diffusion. A prograde torque is 
transmitted from low to high latitudes along these field lines. The slow polar and fast equatorial 
rotation are indicated by the darker shadings of the outer sphere. The dashed lines indicate the 
latitudes at which the rotation of the convection zone matches that of the interior. 



schematically in figure 1(b) [ in which the linkage is via a time-averaged field whose lines 
thread the tachocline, forming the superficial part of a global- scale interior dipole sta- 



bilized by a deep toroidal field (e.g. lBraithwaite &: Spruitil2004[ ). Such an interior dipole 
has a diffusive lifetime somewhat greater than the Sun's lifetime of around 4.5 x lO^yr. 
The dipole imposes an Alfvenic "Ferraro constraint" on the interior. It is this constraint 
that hel ps to enforce the interior's solid rotation (e.g. lFerrarolllQSTtlMestel fc Weisslll987 : 
[MacGr cgor fc Charbonneaulll999l ). 



The field lines shown in figure 1(b) emerge from the interior (light-grey sphere) near 



the north pole and, after threading their way through the tachocline, re-enter the inte- 
rior near the south pole. They return northward through an interior "apple-core" region, 
not shown, surrounding the rotation axis. It is crucial that the field lines emerging from 
the interior bend over toward the horizontal as they enter the tachocline. They must be 
prevented from exte nding upward through the polar cap, as occurs when magnetic diffu- 
sion dominates (e.g. lBraithwaite fc Spruitll2004l : iBrun fc Zahnll2006l ). The curved shapes 



of the field lines in figure 1(b) are evidently such as to transport angular momentum 



from low to high latitudes by means of persistent Alfvenic torques, exactly as required 
to prevent the tachocline's MMCs from burrowing into the interior and thickening the 
tachocline. 



The time averaging envisaged in figure 1(b) conceals a plethora of fast processes, includ 



ing the 22-year dynamo cycle, convective overshoot, and other turbulent processes arising 
from various instabilities in the tachocline. All these are fast relative to the timescales on 
which the mean structure of the tachocline is maintained, ~ lO^yr or more. We presume 
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that the fast processes have two important consequences. The first is to produce a tur- 
bulent magnetic diffusivity that stops the field lines being wound up arbitrarily tightly 
by the shear in the tachocline, keeping the curved shapes shownjj 

The second important consequence is that, away from the poles, the field lines are held 
down, and held approximately horizontal, by turbulent "magnetic flux pumping" from the 
convective overshoot layer. The effectiveness of such fiux pumping can be strongly argued 
from several lines of evidence, including three-dimensional direct numerical simulations, 
with varying emphasis on the role of turbulent anisotropy and of vertical gradients of 
density and turbule nt intensity (e.g. Tobias et aLlbOOll : iKitchatinov fc RiidigeiJuOOd & 



refs.) (see also §3 of lWeiss et aLll2004l for a historical review). 

Near the poles it is less clear that magnetic fiux pumping will be effective in confining 
the field. At least its effectiveness for near-vertical magnetic fields has not, to our knowl- 
edge, been convincingly demonstrated. However, as argued for instance in GM98, there 
are good reasons in any case (iJH below) to expect the tachocline's MMC near the poles 
to take the form of weak but persistent downwelling. This suggests that the field can, 
in any case, be confined in the polar caps through an advection-diffusion balance, the 
kind of balance argued for heuristically in GM98. The purpose of this paper is to show in 
detail, by solving an appropriate set of nonlinear magnetohydrodynamic equations, that 
such polar confinement by downwelling is indeed possible in a physically realistic model, 
applicable to the Sun both today and early in its main-sequence lifetime. 

A large family of axisymmetric nonlinear solutions showing polar confinement has been 
obtained using two different techniques. The first technique is semi-analytical in a sense 
to be explained, and the second is numerical on a 2-dimensional grid. The solutions 
are to be regarded as candidate solutions for possible fiows in the real Sun, all showing 
confinement in the sense that the total magnetic field strength |B| dies off exponentially 
with altitude, thanks to downward advection acting against upward diffusion. In this 
sense the poloidal and toroidal field components are both well confined. We call these 
flows "conflnement layers" . They are not to be confused with the tachocline itself. Rather, 
they occupy relatively thin regions at the bottom of the tachocline and are much more 
weakly sheared, with relatively long timescales ^^ lO^yr. 

The detailed dynamics involves not only magnetic advection, stretching, twisting and 
diffusion but also a near-perfect balance between Lorentz, Coriolis, pressure-gradient 
and buoyancy forces (§Sj3fF.). Thus the confinement-layer flows are magnetostrophic, like 
cert ain flows that have be en studied in connection with models of the Earth's liquid core 



(e.g. iKleeorin et aLlll997l & refs.), though different in most other respects. For instance 



the latter flows are viscous but unstratified: buoyancy forces and thermal diffusion are 
absent. In the confinement-layer fiows studied here, by contrast, viscosity turns out to 
be wholly unimportant while buoyancy and thermal diffusion are crucial, along with 
magnetic diffusion. 

Figure [5] gives a preview of a typical confinement-layer fiow, seen in vertical cross- 
section. It shows the poloidal velocity and magnetic field components from a numerical 
solution. The emerging magnetic field lines are bent over within the confinement layer. 



as required to fit into the global picture sketched in figure 1(b) The magnetic field B 
has a toroidal component, not shown in the figure, imparting spiral shapes to the three- 
dimensional field lines and providing the prograde Alfvenic torque demanded by the 
global picture, in balance with a retrograde Coriolis torque on the equatorward fiow. 

f Of course the persistent angular momentum transport from the curved B lines could be 
supplemented by equally persistent contributions coming directly from MHD-turbulent stresses 
(e.g. ,Spruit„2002. : .Oilman fc Callv„2007. : Jarfrcv fc Menou.,2007. ). 
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Figure 2. The magnetic confinement layer near the north pole in a model for today's Sun. 
The field strength |B| falls off exponentially with altitude z. The toroidal components of B 
and velocity u are not shown. The streamlines with arrows show the downwelling responsible 
for the confinement. If the downwelling were switched off, the field near the p ole would diffuse 
and become nearly vertical, as illustrated for instance in iBrun fc ZahnI (|2006l ). Compositional 
stratification is indicated by shading. The plot is from a numerical solution; the corresponding 
semi-analytical solution looks almost identical. The horizontal and vertical axes are colatitudinal 
distance r and altitude z in units of 5, the advection-diffusion scale defined in (|l.ip . With typical 
parameter values, the scale S is of the order of a fraction of a megametre, ~ 10"'^ of the Sun's 
radius. 



The vertical and colatitudinal distances in figure [5] are shown in units of the magnetic 
advection-diffusion scale 

5 = V/U. (1.1) 

say, where U is the magnitude of the downwelling and r/ is the magnetic difFusivity. 
Throughout this paper, we assume that the confinement-layer flow is laminar, and there- 
fore use molecular or microscopic difFusivity values ([J3|). Issues of stability or instability 
lie beyond the scope of this paper but, close to the pole at least, there appears to be a 
strong case for stability, to be argued in a future paper, arising from the smallness of the 
scale 6. Under reasonable assumptions, d is only a fraction of a megametre, far smaller 
than the thickness of the overlying tachocline which latte r, by contrast, is probably un- 
stable a nd indeed turbulent, as already mentioned (e.g. ISpruiti l2002t iGilman fc Cally 
|2007; Pa rfrevfcMenoiJl2n07h . 

In the present-day Sun's helium settling layer, the top of which corresponds to the 
shaded region in figure [2l a downward gradient of helium concentration reinforces the 
stable stratification due to the sub-adiabatic temperature gradient. Because the diffusiv- 
ity of helium through hydrogen, x ~ 0-9 ^ lO^cm^s"^, is much less than the magnetic 
difFusivity rj w 4.1 x lO^cm^s"^, the helium settling layer is nearly impermeable to the 
confinement-layer flow. Helium advection and diffusion are comparable only in the ex- 
tremely thin "helium sublayer" marked in figure [2] In this and other respects, all the 
solutions in the present p aper supersede those described in a first report on this work 
(|Wood fc Mclntvrd [20071 hereafter WM07). For instance, in WM07 we took x to be 
zero, implying a helium sublayer of vanishing thickness. We also took i', the kinematic 
viscosity, to be zero and allowed a finite slip velocity at the top of the helium settling 
layer, assuming that this slip velocity would in reality be resolved into a weak Ekman 
layer. However, the solutions presented here show that, on the contrary, no Ekman layer 
forms. The slip discontinuity is replaced by a smooth velocity profile across the helium 
sublayer and, as will be shown in ij6l the flow stays essentially inviscid. 



6 T. S. Wood and M. E. Mclntyre 

The plan of the paper is as follows. In f}2] we summarize the reasons for expecting 
persistent downwelling over the poles. In fj3]we present the model equations and in f|4] 
the semi-analytical solutions. Those solutions rely on assuming a self-similar horizontal 
structure that is asymptotically valid in the limit of strong stable stratification. The same 
limit was taken in WM07. 

The validity of the strong-stratification limit is assessed in §Sj5]and[ni which take a thor- 
ough look at the dynamical balances and scalings in the confinement layer and helium 
sublayer respectively. Strong stable stratification means that the thermal and composi- 
tional stratification surfaces are "flat" , meaning gravitationally horizontal, to sufficient 
approximation in some region surrounding the poles, which for reasonable parameter 
values can be quite large in horizontal extent, up to tens of degrees of colatitude. Within 
the helium sublayer, the low magnetic Reynolds number and fiat geometry cause the 
momentum balance to take on the character of flow in a porous medium, as fluid pushes 
horizontally past the field lines. As already indicated, true viscous effects are negligible 
everywhere, even in the sublayer. 

Boundary conditions for the numerical solutions are discussed in [jTl The numerical so- 
lutions themselves are presented and discussed in [J8l They provide cross-checks with the 
semi-analytical solutions plus additional insights. In particular, they directly demonstrate 
the fiatness of the stratification surfaces by solving the full equations, for finite strati- 
fication. The solutions allow the stratification surfaces to tilt as they may, but confirm 
that the departures from flatness are indeed small when the stratification is realistically 
strong. In figure [5J for instance, the departures from fiatness are barely visible. 

In [JS]we discuss a subtlety that arises when comparing the semi-analytical and numer- 
ical solutions in the upper part of the flow. The dynamical balances aloft become delicate 
as the Lorentz and Coriolis forces become vanishingly small. The effects of truncation 
error and other small effects thus complicate the comparison. However, this is something 
of an academic point because of our expectation that, in reality, the confinement-layer 
solutions will need to be matched to a turbulent tachocline aloft, a task that remains a 
challenge for the future. 

In SJTO] we show that the presence of the helium settling layer is not crucial to our 
confinement-layer model. The interior field Bj is sufficient by itself to turn the fiow 
equatorward, and the field remains confined in much the same way. That result has 
relevance to the Sun's early main-sequence evolution. It explains for instance how the 
burrowing tendency could have been held in check from the start, allowing the helium 
settling layer to form. In the concluding discussion, ijlli we consider the implications for 
early solar evolution and lithium depletion. 



2. Downwelling in the polar tachocline 

Our polar-confinement scenario relies on the MMC pattern in the stably-stratified po- 
lar tachocline being robustly and persistently downward above the confinement layer, 
after averaging out any fast fluctuations due to waves and turbulence. As recognized in 
GM98, helioseismology provides a compelling reason to expect polar downwelling rather 
than upwelling, at least in today's tachocline. A further reason is that a downward MMC 
over the pole is a robust consequence of the gyroscopic pumping already mentioned, 
which, in the absence of the interior fleld Bj, would mediate the downward spread- 
ing or burrowing of the convection zone's slow polar rotation. The distinctio n between 
yroscopically-pumped MMCs and MMCs driven in other ways is reviewed in iMcIntvrd 
2007^ §§8.1-8.2), confirming also the robustness of the burrowing tendency itself, de- 
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spite recent controversy. The upshot is that we expect persistent polar downwelhng to be 
present not only in today's Sun, but also throughout the Sun's main-sequence lifetime. 

The argument from helioseismology is as follows. As is well known, the pressure, den- 
sity and angular velocity fields, averaged with respect to time and longitude, satisfy 
hydrostatic and cyclostrophic balance to excellent approximation. Departures from such 
balance must take the form of fast oscillations such as p- modes and g- modes, or turbulent 
fluctuations. From the curl of the momentum equation, taking its azimuthal component, 
we may show in the standard way that balance implies the so-called "thermal-wind re- 
lation". In cylindrical polar coordinates (z,r, 0) centered on the axis of rotation, with 
the axial coordinate z directed vertically upward at the north pole, the thermal-wind 
relation can be expressed as 

rp^^}^ = (y^py,'^p).e^ (2.1) 



where J7 is the absolute angular velocity of the Sun's differential rotation, p is the density, 
and p is the total pressure. The unit vector e^ is directed azimuthally, while Vp, being 
dominated by its hydrostatic part, is very close to being vertically downward. On the 
assumption that the observed negative sign of d\n\'^/dz persists into the region near the 
pole invisible to helioseismology — Occam's razor makes this a reasonable assumption 
— we must have a minimum in p, and hence a maximum in temperature T, on each 
isobaric surface at the pole. 

The stably-stratified radiative envelope is a thermally relaxing system. Local temper- 
ature anomalies, defined as departures of T from local radiative equilibrium, will tend to 
relax back toward zero. To hold T above radiative equilibrium near the pole, there has to 
be persistent adiabatic compression by downwelhng, with compensating upwelling and 
negative T anomalies in lower latitudes. 

The strength U of the polar downwelhng is difficult to estimate precisely. Among 
other things it depends on the tachocline thickness, which is not well constrained by 
helioseismology. The thickness scale enters both via equation (|2.1I) and, more sensitively, 
via the rate of diffusive thermal relaxation within the tachocline. GM98 estimated U ~ 
10~^cms~^, using the rather small tachocline thickness estimate, 13 Mm, derived by 
Elliott fc Goughl ( 19991 ). The value of U thus estimated is inversely proportional to the 



cube of the tachocline thickness, and so a similar estimate using a deeper tachocline 
would yield a much smaller value of U. 



However, GM98 assumed that the bulk of the tachocline is laminar. iMcIntvrd (J200l1) 
considered an alternative scenario in which magnetohydrodynamic turbulent stresses 
within a deeper tachocline dominate the angular-momentum transport from the overly- 
ing convection zone, except near the bottom of the tachocline. The turbulent stresses 
were e stimated by assuming a particular prescription for the turbulence, following iSpruiti 



(|2002l ). The stresses diverge in a thin layer near the bottom of the tachocline, just above 
the confinement layer, where they gyroscopically pump a downwelhng of the order of 
[/ ~ 4 X lO^^cms^^ or greater. 

Fortunately, our confinement-layer solutions can accommodate a wide range of uncer- 
tainty over the value of U. They will show that polar field confinement by downwelhng 
is robust over a range of U values at least as wide as 10~^cms~^ to 10~''cms~^. From 
here on we use GM98's value U ~ 10~^cms^^ for illustrative purposes. 
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3. The model equations 

Consider the magnetic confinement layer near the north pole. As already noted, the 
magnetic advection-diffusion thickness scale 6 = ij/U is to be evaluated with the mi- 
croscopic magnetic dif fusivity y, who se value in the neighbourhood of the tachocline is 



carefully estimated bv iGoughl (|2007l ) to be 77 w 4.1 x lO^cm^s ^. This gives the value 
6 « 0.4 Mm if C/ « 10-5 cm s-^ 

We work in a frame rotating with the same angular velocity as the interior, fli = 
2.7 X IQ-^s"^, and seek axisymmetric solutions of the Boussinesq MHD equations within 
a domain consisting of a cylindrical volume V surrounding the pole. Cylindrical coor- 
dinates (z, r, (j)) centred on the rotation axis will be used, with corresponding unit vec- 
tors (e^, Gr, e^). The use of cylindrical coordinates will lead to significant mathematical 
simplifications. We may regard them as slightly-distorted spherical coordinates, with r 
representing approximate colatitudinal distance and z locally vertical. The confinement- 
layer flows to be studied are thin-shell flows, with z ^ 5 and r ^ S, and so the coordinate 
distortions should be qualitatively unimportant out to colatitudes as far as 20° or so. 

The Boussinesq framework should itself be highly accurate because typical flow and 
Alfven velocities are tiny fractions of the local sound speed, and because S values of 
the order of a fraction of a megametre are far smaller than the pressure scale height, 
« 60 Mm. Conveniently, Boussinesq dynamics permits us to measure the strength of the 
magnetic held B in terms of the corresponding Alfven speed, with 1 gauss corresponding 



to 0.6cms ^ at a (constant) tachocline density of 0.2gcm ^ (jGoughll2007[ ). 

We impose uniform downwelling of magnitude U aloft and a simple axial dipolar, fully- 
diffused poloidal magnetic field structure beneath, to represent the interior magnetic 
field Bi = {Biz, Bir, Bitj,), on to which the field B in the confinement layer is to be 
matched. This interior dipolar field has B^r/r constant and Biz a linear function of z. It 
is possible to have Bicf, 7^ with B\^jr constant; however, for reasons to be explained at 
the end of the section, the main focus will be on the purely poloidal case B^^ = 0. This 
implies the vanishing of the Alfvenic torque beneath. For the semi-analytical solutions, 
the condition B\^ = is imposed directly. For the numerical solutions a less direct 
procedure is necessary, to be explained in fJS] 

It is convenient to nondimensionalize the equations using b as the lengthscale in the 
horizontal (r) as well as in the vertical (z) direction. Thus the thin-shell nature of the 
fiow is expressed by the dimensionless relation r ^ 1. We take U as the scale for the 
velocity field u, and SjU as the timescale, ~ lO^yr ii U ^ lO^^cms^^. This is the 
advection timescale for the fiow through the confinement layer and, by construction, is 
also the timescale on which B diffuses across the confinement layer. Since this timescale 
far exceeds the typical turnover time of the turbulent eddies in the overlying layers, we 
may neglect any fiuctuations in the downwelling aloft. We therefore take \J to be steady 
as well as uniform, representing the time-averaged downwelling that is gyroscopically 
pumped by turbulence in the overlying layers, whether those layers consist of the convec- 
tion zone or the tachocline or both. Fluctuations in V may well be present, but should 
not greatly infiuence the structure of the confinement layer on timescales ^ 10"'yr or 
more. Even if the downwelling comes in pulses, the cumulative effect will arguably be 
much the same as if it were steady. 

We nondimensionalize the magnetic field B (expressed as Alfven velocity) with respect 
to a different velocity scale (%l\r\)^l'^ k, 0.05cms^^. The significance of this last choice 
will emerge in fj5l It will simplify the scaling relations (I5.5|) - (|5.10[) . We suppose that the 
thermal and compositional stratifications are approximately uniform within the helium 
settling layer, shaded in figure [21 since the confinement layer's MMCs do not penetrate 
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into that region. Writing -d and /t for the fractional (therefore dimensionless) pertur- 
bations of Boussinesq potential temperature and mean molecular weight, we therefore 
impose that the corresponding buoyancy frequencies are exactly constant at the bottom 
of the domain. That is, we impose, with z now the dimensionless vertical coordinate, 



d-d 
dz 

dfi 
dz 



NI6 



= const., 



bottom 



bottom 



9 



const. 



(3.1) 
(3.2) 



the dimensional buoyancy frequencies N{) and Nn bein g constant by definition. For to- 
day's Sun we have N^ w 0.8 x 10~^s~^ (IGoughl 120071 and Gough 2010, personal com - 
munication), and N^ w 0.5 x lO^'^s^^ (e.g. IChristensen-Dalsgaard fc Thompson! l2007l ltl 



corresponding to a total buoyancy frequency N = {N^ 



^2)1/2 



0.94 X lO-^s-i rep- 



resentative of the stratification just inside the helium settling layer. In place of i? and p, 
it proves convenient to define rescaled quantities d and fi, also dimensionless, by 



^, = k 



9 



fi = II; 



(3.3) 
(3.4) 



so that the dimensionless stratifications inside the helium settling layer become simply 
d-d/dz = 1 and d^/dz = — 1. 

Finally, we nondimensionalize the pressure anomaly p' by 2flirip where p « 0.2gcm^'^, 
the constant Boussinesq density, and thus arrive at the following dimensionless equations. 



Ro— + e^ X u = - Vp' + a.Q-de^ - a^^e^ 

+ (V X B) X B -h EkV^u 
= V-u 

^- = Vx (uxB) + v2B 
at 

= V-B 

^ = ^V^^ 
Di 77 

Dm X^2 
Bt ri ^ 



(3.5) 
(3.6) 
(3.7) 
(3.8) 
(3.9) 

(3.10) 



where D/Dt — d/dt + u • V, the material derivative, and where e^ is the vertical unit 
vector. The thermal and compositional diffusivities k and x take numerical val ues k 



1.4x 10^ cm^s ^ and X ~ 0.9x 10^ cm^s ^ in the neighbourhood of the tachocline (JGough 



20071) . In (13.101) we have neglected gravitational settling, an excellent approximation in 
virtue of the short timescale ~ lO^yr of the confinement-layer dynamics relative to the 



f Estimates 



of 
-3^ 



A'^ vary (e.g. IChristensen-D alsgaard et al\ Il993f ). The value 
Np, ft; 0.5 x 10~'^s~^ was comput e d in Mclntyrc (2003, §8.5) from information given in 
IChristensen-Dalsgaard fc ThompsonI (|2007l ). However, the results presented in this paper are 

see iJSl 



not critically dependent on the value of A^^ 



(3.11) 
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Sun's lifetime, > lO^yr. We have defined four dimensionless constants 

where v is the kinematic viscosity, the sum of molecular and radiative contributions, 
« 2.7 X lO^cm^s^^ in the neighbourhood of the tachocline (Gough 2007). The Rossby 
and Ekman numbers, Ro and Ek, quantify how far magnetic flux and fluid momentum 
diffuse across the confinement layer during one solar rotation. For V ^ lO^^cms^^ the 
Rossby number is tiny, Ro ^ 0.5 x 10^^, and the Ekman number is smaller still because 
V jj] K, 0.7 X 10^^. To excellent approximation, therefore, the fiows under consideration 
will be magnetostrophic. That is, in p.Sp the Coriolis force will be balanced against the 
combined pressure-gradient, buoyancy, and Lorentz forces: 

e^ X u = -Vy + oi^-dez - a^/^e^ + (V x B) x B , (3.13) 

and this will be verified independently from the numerical solutions, from which magne- 
tostrophic balance emerges rather than being imposed. 

We are concerned here only with axisymmetric steady states. Then the azimuthal 
components of (13.13^ and its curl are respectively 

M, = ^B • V(rB^) , (3.14) 

r 



a 



dz dr dr r d 



z 

where r is the dimensionless perpendicular distance from the rotation axis, and where 
suffixes z, r, cf) denote vector components. 

Equation p.l4p represents the local torque balance about the rotation axis, after mul- 
tiplication by r. It describes how the retrograde Coriolis torque from the equatorward 
flow is balanced by the prograde Lorentz torque from the confined magnetic fieldljj Equa- 
tion (|3.15p represents thermal-wind balance generalized to include compositional gradi- 
ents and the Lorentz force-curl. In the upper part of figure [2l where the magnetic field 
and compositional stratification are both negligible, this balance becomes the standard 
thermal-wind balance „ „ „ 

^=«.^> (3-16) 

which is the Boussinesq, low-Rossby-number limit of (j2.ip in dimensionless units. 

The overall torque balance for the confinement layer can be expressed by integrating 
r times p.l4p over the volume V of the cylindrical domain, then using the divergence 
theorem and the fact that V • u = 0. The result is 

ir^u-dS == ( rS^B-dS (3.17) 

dV JdV 

where dS is the vector area element directed outward. In the corresponding dimensional 

f A referee reminds us that this torq ue balance is related to the well-known Taylor constraint 
for magnetostrophic flow (jTavlorll 19631 ). If (|3.14p is integrated vertically between two hypothet- 
ical impermeable boundaries then, because V ■ u = 0, the integral on the left and hence that on 
the right must vanish. In the conflnement-layer problem, however, this constraint is broken by 
the presence of downwelling aloft. 
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equation, ^r^ is replaced by ftir^, the absolute angular momentuni per unit mass after 
neglecting contributions 0(Ro). 

The right-hand side of (|3.17p represents the total Alfvenic torque exerted on the con- 
finement layer. The left-hand side represents the net rate of absolute angular momentum 
export by the flow coining in through the top and out through the periphery, as in 
figure El For any such velocity field the left-hand side is positive. Therefore the Alfvenic 
torque must also be positive, i.e. prograde. 

The need for the field Bi is now apparent. Without it, the balance described by p.l7p 
would be impossible. Instead, the fluid within V, and outside it as well, would begin to be 
spun down by the flow. In the Sun, this would cause the slow rotation of the high-latitude 
convection zone to spread down into the radiative envelope, in the Haynes-Spiegel-Zahn 
burrowing process noted in fJTl We also note that the torque balance described by ( 3.17P 



is very d ifferent from th e torqu e balance described bv lRiidiger fc Kitchatinovl (|l997l ) and 



iKitchatinov fc Riidigeij ( 20061 ) in which there is no Coriolis effect, the Lorentz torque 



being balanced instead by a viscous torque. In fact for realistic solar parameter values it 
will be seen that, as already indicated, viscous torques are entirely negligible — perhaps 
counterintuitively for shear layers as thin as the confinement layer and the helium sub- 
layer. 

The balance expressed by p.l7p must apply also throughout the interior apple-core 
region that surrounds the rotation axis, magnetically linking the confinement layer at the 
north pole to that at the south. The mu-choke, both in the Sun's helium-rich inner core 
and in the helium settling layer, if present, is enough to suppress MMCs in the apple-core 
region and make the left-hand side of p.l7p negligible, when the volume of integration 
V is taken as the interior apple-core region. Any Alfvenic torque exerted on the bottom 
of one confinement layer must therefore be balanced by an equal and opposite Alfvenic 
torque exerted on the bottom of the other confinement layer. If we assume that the two 
confinement layers are mirror-symmetric about the equatorial plane, then it follows that 
the torque, and therefore Bj^, must vanish on the bottom of each confinement layer. The 
prograde Alfvenic torque on each confinement layer must therefore come wholly from its 
sideways connection to lower latitudes via the tachocline, which is consistent with the 



global picture suggested in figure 1(b) 



4. The semi-analytical solutions 

Equation p.lSp describes how the Coriolis and Lorentz force-curls act to tilt the strat- 
ification surfaces within the confinement layer. Now if the stable stratification is suffi- 
ciently strong, then the tilting will be only slight. The formal asymptotic limit to describe 
this is a^ , a^ — ;■ cx) with u and B finite so as to preserve the steady-state torque balance 
(|3.17p and a steady-state balance in the induction equation (13. 7p . In that limit, (|3.15l) im- 
plies that dd/dr -^ and d^/dr — )► 0, with both sides of (|3.15p finite. The stratification 
surfaces become perfectly flat: 'd — >■ 'd{z) and, in the helium sublayer, /x — )► ^l{z). 

In the limit of perfect flatness thus enforced by p.lSp . the remaining equations can 
be reduced to a system of coupled ordinary differential equations. With § — ^{z) and 
H — ^[z), (13. 9p and (|3.10p imply for steady flow that 

-—In— = -—In— (4 1) 

rj Az Az Tj Az Az 

Therefore Uz is also a function of z alone, Uz = Uz{z). From (j3.6p it then follows that Uj. 
is r times a function of z alone, on the assumption of regularity at the pole r = 0. We 
say that the poloidal velocity field is "horizontally self-similar" . The induction equation 
(|3.7p then permits a steady poloidal magnetic field that is horizontally self-similar in the 
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same sense. The resulting equations are 



rUr^B- V{rB^) (4.2) 

= 1^ + ^ (4.3) 
r or az 

..V^^-rB.V^+fv^-l'ji?^ (4.4) 
r r \ r'^ J 

0^^^^ + ^. (4.6) 



r 



dr dz 



The azimuthal component of p.Sp is replaced by p.l4p and multiplied by r to give 
expressing the balance between Coriolis and Lorentz torques as before. Equations (|4.3p - 
(|4.6p correspond to p.6p - p.8p : equations p.9p and (I3.10p have no further role, beyond 
their connection to the downwelling expressed by (j4.ll) . Equation (|4.4p permits more 
general toroidal magnetic and differential-rotation fields than were considered in WM07. 

For large but finite a^ and a^, departures from perfect flatness arise as small correc- 
tions. To describe these corrections it is necessary to bring back equations (|3.9|) . p.lOp 
and (|3.15p : the order of magnitude of the corrections is analysed in §iJSHni The whole 
picture will be independently checked by the numerical solutions, which automatically 
contain the departures from flatness governed by (I3.9p . p.lOp and (I3.15P since the full 
set of equations is used, with finite values of a^ and a^t, for instance to produce the 
almost-fiat numerical solution plotted in figure [2] 

Returning now to the limit of perfect flatness, we focus on equations (I4.2p - (|4.6p . These 
equations admit a family of solutions in which the function Uz{z) is arbitrary except for 
certain restrictions on its asymptotic behaviour as z — > ±cx). In particular, we require 
Uz{z) -^ —1 as z -^■ -l-oo and Uz{z) — > as z ^- — oo. These statements will shortly be 
made more precise. As discussed in §7, the arbitrariness in Uz{z) is needed to permit 
matching to the confinement layer's surroundings. 

Solutions can be most conveniently constructed by taking advantage of the arbitrari- 
ness to specify a suitable Uz{z) at the start. With Uz{z) specified, we can then find Bz{z) 
numerically by solving the vertical component (|4.5I) of the induction equation as a lin- 
ear ordinary differential equation, assuming that the time-averaged field Bz vanishes far 
above the confinement layer (z -^ -|-cx)) and that it matches on to the imposed interior 
dipolar magnetic field structure beneath {z — s- — oo). The interior dipole has the same 
horizontally self-similar structure as the confinement layer, with components satisfying 
Bir/f = constant, and B-^z a linear function of z consistent with (14. 6p . Even though the 
balance in (|4.5p is not simple advective-diffusive, we find that Bz decays upward like 
exp(— z). 

The radial components of u and B can be found directly from their vertical compo- 
nents, by using ()4.3I1 and ()4.6p and assuming regularity at the pole r = 0: 

"--^^' («, 

---5^^ (.8, 

So once we have Bz{z) we can calculate Br from (|4.8I) . and then the toroidal field B^ 
from ()4.2p by using ()4.7p and taking advantage of the hyperbolic character of the operator 
B • V. By calculating B^ in this way, we ensure that the Lorentz torque balances the 
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Coriolis torque along each magnetic field line. Requiring that B^{r, z) — > for all r as 
z -^ —oo (recall the end of ^ leads to the following, unique solution of (|4.2p : 



B^ ^ bJ ^dz. (4.9) 

For any Ur profile that decays exponentially as z ^- — oo, this solution for B^, and with 
it the Maxwell stress and Alfvenic torque, will also decay exponentially as z — > — cxd. The 
expression (14. 9p then shows that B^ has the same horizontally self-similar functional 
form as Ur and Br, namely r times a function of z alone. 

To ensure that B^ decays aloft, as z — > +cxd, it is sufficient to assume that 

Ur — 0(exp(— 7z)) as z ^)- +00 (4-10) 

for constant 7 > 1, implying that Uz(z) ~ — l-|-0(exp(— 7z)). The assumption that 7 > 1 
ensures that we get solutions with qualitatively reasonable behaviour aloft. Otherwise we 
leave the value of 7 arbitrary. Once again this is an arbitrariness whose resolution will 
depend on matching to the confinement layer's surroundings, in this case to conditions 
aloft. As already indicated, the conditions aloft probably involve turbulent flow, for which 
we do not yet have quantitative models. So here we restrict ourselves to surveying the 
possible range of behaviours for 7 > 1 . 

The three cases 7 > 2, 7 = 2, and 2 > 7 > 1 need separate consideration. When 
7 > 2, the only case considered in WM07, the integral in (|4.9|) converges to a constant 
plus 0(exp(— (7 — 2)z)) as z — > +cx). That in turn means that B^ decays upward like 
exp(— z). When 7 = 2, the integral in ()4.9p asymptotes to a linear function of z, and B^ 
decays upward like zexp(— z). When 2 > 7 > 1, the integral in (14. 9p increases upward 
like exp((2 — 7)z), but B^ still decays upward, like exp(— (7 — l)z). 

In all these cases it is clear from ()4.9|) that B^f, does not have exactly the same z- 
dependence aloft as does Bz- This is contrary to what might might have been expected 
from a naive appeal to advective-diffusive balance, with advection by constant down- 
welling Uz = —1. Having the same z-dependences would make the right-hand side of 
(|4.2|) vanish. Hence it is the more or less subtle departures from advective-diffusive bal- 
ance, including the contribution to B^ from the twisting of field lines by the differential 
rotation u^, that allow the right-hand side of (|4.2[) not to vanish and thereby to provide 
a Lorentz torque to support the flow Ur at all altitudes. 

The U0 fleld that does the twisting can be calculated next, from (|4.4I) and the condition 
that the interior rotates solidly, w^ — > as z — )■ —00. Again this calculation depends on 
the hyperbolic character of B • V. When B^ is given by (|4.9p we have, uniquely, 

dB,_d^\dz_ 

" dz dz^ J Bz ^ ' 

showing that m^ is also r times a function of z alone. That is, the differential rotation is 
what astrophysicists call "shellular solid rotation". In the three cases 7 > 2, 7 = 2, and 
2 > 7 > 1, the behaviours of u^ as z — > -l-oo are respectively u^ ^ constant, u^ ^ ±z, 
and U0 ~ ±exp((2 — 7)z). 

Cases with negative shear aloft, du^/dz < — especially the last case, with exponenti- 
ally-increasing negative shear aloft — are suggestive of a possible way to match upwards 
to the observed, much stronger negative shear in the bulk of the tachocline. By using 
([49]) to eliminate B^ from d/dz of (|4lT|) . then (|4?5l) to eliminate d^B^/dz^ and ([iJ]) to 
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Figure 3. Vertical profiles from a semi-analytical solution of the confinement-layer equations 
in the strong-stratification limit q^, q^ — >■ oo. The downwelling profile Uz{z), solid curve on the 
left, was chosen to match the downwelling profile from the numerical solution shown in figure (2] 
For numerical reasons, small adjustments were made to this profile in the "slippery" upper 
region z > 1.55; see §3IH9l In (|4.10p the decay constant 7 — 2.24, and the u^ profile therefore 
approaches a constant like exp(— 0.24z). The parameter A, (|5.4p below, takes the value A ~ 3.5. 



eliminate duz/dz, we find 



dz 



UzUr 



1 dUr 



2Ur 



S? 



/ in"'' 

\1~^)~fy? ^^ z — > 00 . 



(4.12) 



The asymptotic behaviour on the right comes from the first two terms in the exact 
expression. The third term involving the integral is smaller by a factor 0(exp(— 72;)). 
Asymptotically, therefore, the sign of the shear du^/dz aloft is the same as the sign 
of Ur aloft. We can therefore find solutions that match on to the strong negative ta- 
chocline shear provided that there is an exponentially weak poleward mass flux above 
the confinement layer. However, a more precise description of such matching must await 
future work, for reasons already mentioned. We do not yet have quantitative models of 
the precise conditions aloft, which are likely to be affected by sma ll-scale MHD turbulence 
(e.g. ISpruitl[200llGilman fc Callvl[2007t IParfrev fc Menoul[2007l & refs.). Aspects of this 
are touched on again in ^ the main point for present purposes being that the structure 
aloft is sensitive to conditions aloft whereas, as is clear from (j4.9|) and (|4.11|) . the rest of 
the confinement layer is insensitive to conditions aloft, as we have verified by varying 7. 

Purely for illustration we show one of the solutions in figure [31 somewhat arbitrarily 
choosing 7 = 2.24. In this case the interior field Bj is taken such that B^/r = 1. The 
downwelling profile Uz{z) was adapted from the numerical solution shown in figure [21 in 
the manner described in S|9l 

Some three-dimensional streamlines and magnetic field lines corresponding to the so- 
lution in figure [3] are plotted in figure [H visualizing how the prograde Lorentz torque on 
the right of ()4.2|) . associated with field- line curvature, balances the retrograde Coriolis 
torque on the left of (|4.2|) and satisfies the overall torque balance p.l7p . 
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Figure 4. Streamlines (left) and magnetic field lines (right) from the semi-analytical solution 
whose vertical profiles are displayed in figure |3l The peripheral shading locates the helium sett- 
ling layer and helium sublayer. 



5. Confinement-layer scalings 

We are interested in the departures from perfect flatness associated with large but 
finite stratification. We may regard those departures as the 0(e) correction terms in an 
asymptotic expansion whose leading, 0(1) term is a semi-analytical solution, where the 
small parameter e is inversely proportional to the stratification and where the limit e — >■ 
is taken within a cylindrical domain of fixed dimensionless size r = r^. 

ft will be found that the departures from flatness behave like 0{er^) near the pole. So 
instead of saying that the departures are small within the fixed domain r ^ r^ we may 
say, perhaps more usefully, that the semi-analytical solutions are valid as leading-order 
approximations as long as we are well within some dimensionless colatitudinal distance 
of the pole, r^ say, that is large in comparison with unity. For a fair range of parameter 
values the dimensional counterpart of r^ turns out to be quite large numerically, of 
the order of hundreds of megametres, or tens of degrees of colatitude. Of course for the 
solutions to apply we must also be within a region of approximately uniform downwelling. 
In this section we use scaling arguments to arrive at an appropriate definition of r^ for 
the bulk of the confinement layer, where the stratification is purely thermal. The next 
section presents the corresponding analysis for the helium sublayer. 

Consider, then, the scaling regime in the bulk of the confinement layer. Because the 
photon mean free path makes the thermal diffusivity k relatively large, with k/tj ^ 3x 10'*, 
the confinement-layer flow only weakly perturbs the background thermal stratification. 
Consistently with (13. 1|) and p.3|) we define a dimensionless thermal anomaly i}' such 
that -d = z + •&' + constant, and such that f?' — ?> beneath the confinement layer at 
the pole. Then, at the pole, and by implication sufficiently close to it, the leading-order 
balance in the dimensionless steady-state thermal equation (|3.9|) involves only the vertical 
component Uz of u, 

uz = -w'^r . 

V 
On the right-hand side we have V^i9' '^ d^d' /dz^ 
d/dz ^ 1, and since ■(?'—>■ beneath. Hence with \u 



(5.1) 

•d' since in the confinement layer 
-^ 1 we have 



i3' - ry/K w 3 X 10"^ < 1 



(5.2) 



In the bulk of the confinement layer we may estimate the departure from fiatness, equiv- 
alently the r-dependence of d' , from (13.151) with a^ neglected and a^ considered large. 
Since the leading-order solution is a semi- analytical solution, the remaining terms in 
(|3.15p all take the form r times a function of z alone, to leading order, from which we 
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may deduce that dd' /dr is r times a function of z alone and hence that 

^' = a^ Jr b^r\ (5.3) 

to leading order, where a^ and b^ are functions of z alone and where a^ has the small 
magnitude given by (I5.2p . a^ ~ tj/k. With d/dz ~ 1 the condition for (|5.3|) to be 
compatible with uniform downwelling, Uz = Uz{z) in (j5.1|) . can now be seen to be r <C 
lan/b^l^/"^. We may therefore take r^ ~ \ai)/b^\^^'^ or, alternatively, e ^ \b§r^/a^\. The 
magnitude of b^ is governed by the remaining terms in p.lSp . 

Perhaps surprisingly, the magnitudes of b^ and r^ can be simply related to a single 
magnitude, that of the vertical component Bz of the magnetic field. This is because the 
velocity and magnetic fields of the leading-order, semi-analytical solutions have compo- 
nents that are all simply related to Bz, thanks to the horizontally self-similar structure. 
Let B be the dimensional magnitude of Bz at the bottom of the confinement layer, in 
units of Alfven speed. Then the corresponding dimensionless magnitude of Bz is A^'^, 
where 

p2 

A = ;-^ , (5.4) 

an Elsasser number based on B. We assume that A^/^ characterizes the order of magnitude 
of Bz in the bulk of the confinement layer, and that d/dz '--^ 1 as before|jJ Then the 
leading-order, semi-analytical solutions satisfy the dimensionless scaling relations 

(5.5) 
(5.6) 
(5.7) 
(5.8) 
(5.9) 
(5.10) 

and these relations are expected to apply at least for r <^ r^, and probably also for 
r < r^ as a guide to orders of magnitude. They can be derived alternatively from scaling 
arguments applied directly to p.l4p along with the steady-state versions of (|3.6p -- p.8p . 
We note that, since the horizontal components of B increase with r, the total magnetic 
field strength |B| can greatly exceed B for r S> 1. 

For increasing values of A the field lines become stifFer, so that both they and the 
velocity streamlines spiral less tightly, as noted in WM07. The magnitude of A is not 
well constrained by observations, depending as it does on the magnitude of the interior 
field at the top of the radiative envelope. Fortunately, however, it will be found that 
the confinement-layer regime described here can accommodate a considerable range of A 
values. 

Now in p.l5p . with the term in a^ neglected, the term in a^ cannot exceed the largest 
of the other terms in order of magnitude. From (15. 3p we have d'd' /dr ~ b^r so that 
the term in a^ has magnitude ^ a^b^r, at least for r ^ r^- The remaining terms in 
p.l5|) have magnitudes either ~ rA~^ (the terms in u^ and S?) or ~ rA (the last 
term). Those magnitudes follow from the horizontally self-similar structure of the semi- 
analytical solutions, along with d/dz ~ 1 and the magnitudes (|5.5p - (|5.10|) . We may 

t These scaling assumptions are consistent with the example solution shown in figure |3l for 
which we imposed Br/r = 1 below the confinement layer, implying that B^ is of order unity 
within and just below the confinement layer. Taking 2 = as the bottom of the confinement 
layer in figure|3l we read off Bz = A^" x 1.9. Consistently, the numerical output gives A ~ 3.5. 



Uz '- 


- 1 


Ur '■ 


■^ r 


U0 - 


^rA^i 


Bz - 


^Ai/2 


Brr 


^ rAi/2 


B^r 


^ rA-1/2 



Confinement of the Sun's interior magnetic field 17 

therefore define the typical magnitude of b^ to be a^ max(A,A~^). Correspondingly, 
with a,3 '^ rj/n we may define r|, the typical magnitude of a^/b^, to be 

-2 - '''''^.min(A,A-i) = M^ min(A, A"!) , (5.11) 



K 



2n^K. 



since a^ = A^J(5-^/2r2i7y. For realistic A'^^ w O.SxlO^'^s^^, for downwelling C/ ~ lO^^cms^^, 
and for A ^ 1, we have 5 ~ 0.4Mm and r^ ■^ 4 x 10"^, corresponding to quite a large 
dimensional colatitudinal distance r^6 ^ 1600Mm. |jj The validity of the foregoing scale 
analysis, and that of the next section, will be independently checked by the numerical 
solutions in SJS] 

Of course our cylindrical model with its assumption of uniform downwelling will itself 
cease to apply, almost certainly, well inside such large distances from the pole. One rea- 
son for this limitation is that the cylindrical coordinates, regarded as distorted spherical 
coordinates, become increasingly inaccurate at large distances from the pole. However, 
that is not the most serious limitation. Using true spherical coordinates would not qual- 
itatively alter the dynamics of the confinement layer. It is the assumption of uniform 
downwelling that places the more serious limitation on the range of applicability of our 
model. At some colatitude the downwelling from the real tachocline must give way to 
upwelling, as required by mass conservation. The confinement- layer regime cannot then 
apply even qualitatively. Instead, the interior magnetic field lines are free to advect and 
diffuse upward until they encounter the magnetic flux pumping associated with the con- 



vective overshoot layer, as assumed in figure 1(b) This has wider implications to be 
discussed in ijlli including implications for lithium burning. 

The range of interior field strengths accommodated by the confinement-layer regime 
is determined by (jS.lip . For uniform downwelling U, the condition for the regime to 
apply quantitatively within, say, 10° colatitude or 90 Mm of the poles is r^S ^ 90 Mm. 
We assume qualitative applicability for r^S > 90 Mm. We can use (I5.11|) together with 
realistic N^ and diffusivity values to write this last condition as 



c(A, A-i) < 3 X 10^ ( i^ y] 



(5.12) 



So for U ^ lO^^cms^^ we expect the regime to apply qualitatively over a range of more 
than four decimal orders of magnitude in A. The corresponding range of field strengths, 
being proportional to A^' ^, covers more than two decimal orders of magnitude. To relate 
this to global-scale interior field strengths, we estimate \Br\ from (|5.9|) at a nominal 
30° colatitude, with S ~ 0.4 Mm so that r « 650. Converting to dimensional units, we 
see that the range of field strengths is roughly 3 gauss < \Br\ < 10'^ gauss, near the 
top of the radiative envelope at 30° colatitude. It is noteworthy that these values lie 
substantially above the threshold, more like 10~^ gauss, for the field strength required to 
enforce the Ferraro constraint in the inte rior over the Sun's lifetime (e.g. Mestel fc Weisj 



1987t ICharbonneau fc MacGregoiJll993f ). For smaller values of U, a still wider range of 



interior field strengths becomes possible. 

f Even when the tilting is significant, such as to become incompatible with (|5.ip . the slopes of 
the thermal stratification surfaces are still geometrically small — far smaller than the geometrical 
aspect ratio r^ . Indeed, even on a global scale we expect the stratification surfaces to depart 
from the horizontal by only "a very tiny fraction of a megametre" from pole to equator (jMcIntyrd 
l2007l . end of §8.5), based on observational constraints on shear in the tachocline. Here of course 
"horizontal" means tangential to the heliopotentials, i.e. to the sum of the centrifugal and 
gravitational potentials. 
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6. The helium sublayer 

In this section we show that the flow through the hehum sublayer has the character 
of flow through a porous medium. We also show that, over a large range of N^ and a^ 
values, the tilting of the compositional stratification surfaces in the sublayer is even less 
significant than that of the thermal stratification surfaces in the confinement layer. 

The helium sublayer marks the transition between the compositionally well- ventilated 
confinement layer and the nearly impermeable, compositionally stratified helium settling 
layer. Therefore, we expect the dimensional sublayer thickness scale, 6y^ say, to be deter- 
mined by a balance between advection and diffusion of helium. In (I3.10p the advection 
operator u • V scales like the strain rate ^ U/6 just above and within the sublayer, 
because w^ — !> just beneath. The strain rate U/6 must therefore be comparable to the 
helium diffusion rate x/^x' ^i^^e U/6 = 'q/6'^, 

6x - {x/vf'^5 « i<5 (6.1) 

for realistic solar parameters. 

Although the sublayer thickness scale 6x is small in magnitude relative to (5, it never- 

1/2 

theless greatly exceeds the Ekman thickness scale fek = {v / IVl-/) . Specifically, 





X6^ 2n, xS 217i 
77 1/ U ly 


r 


- -Ro-i > 1, 



(6.2) 

because x/'^ ~ 0-3 while Ro^'^ ^ 1, typically by many decimal orders of magnitude; recall 
(|3.12p . The relations (I6.ip and (|6.2[) suggest that the dynamics of the helium sublayer 
should be well described by the asymptotic regime 

fek < <5x < (5 . (6.3) 

We assume (|6.3p throughout this section. 

Under (|6.3p the magnetic diffusion rate v/^x ^^ ^^"^ sublayer greatly exceeds the helium 
diffusion rate x/ ^\i by a factor r\/x ^ (^/^x)^- "^^^ flow within the sublayer can therefore 
induce only a small perturbation B — Bj = B', say, to the interior fleld Bi. In figures [H 
and m the field lines are hardly deflected as they cross the sublayer. We may therefore 
analyse the sublayer as a perturbation to the state with u = and B = Bi, where B; 
has the simple dipolar structure already assumed, with components satisfying B\^ = 0, 
Bir/r = constant, and Biz a linear function of z consistent with V-Bi =0. 

Any such Bi has V x Bi = and Lorentz force (V x Bi) x Bi = 0. Using this we show 

in Appendix lAl that, in the asymptotic regime given by ()6.3|) . the steady-state induction 

equation becomes simply „ „2 

- Bi,— u+-^B', (6.4) 

oz oz'^ 

in the dimensionless variables of fJSl The momentum balance (|3.13p becomes 

d_ 
' dz 

Here the thermal-buoyancy term a^iJe, has been absorbed into a modifled pressure 
gradient Vp, which also incorporates a gradient contribution to the Lorentz force; see 
Appendix El below (IA3p . Because 6x -^6 and because, as verifled shortly, the sublayer 
will prove to be sufficiently flat, we may take B^z to be constant throughout the sublayer. 
It is convenient to equate the dimensional value of this constant to B in the definition 



e, X u = -Vp-a^^e, +Bi2^--B'. (6.5) 
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(j5.4p of the Elsasser number A. Then the diniensionless magnitude of Biz in the sublayer 
is precisely A^/^. We can now integrate (|6.4p to give 

1 dB' 

= A^ u + ^^ (6.6) 

oz 

since both B' and u vanish beneath the sublayer. Using (|6.6p we write (16.51) as 

e^ X u = — Vp — a^/ie^ — Au. {^■'^) 

The term — Au has the form of a Darcy or Rayleigh drag, showing that the sublayer 
behaves like a porous medium on the timescale set by the strain flow. The impermeability 
of the helium settling layer together with the sublayer's flatness and thinness act to keep 
the flow nearly horizontal, compelling it to push past, and slightly deflect, the magnetic 
field lines spanning the sublayer at angles steep by comparison with sublayer aspect 
ratios. So the Lorentz force from the deflected field lines resists the sublayer flow in the 



manner of a Darcy friction; similar behaviour occurs in Hartmann layers (e.g. iDebnath 



19731) . although in such cases viscosity also contributes to the balance of forces. Within 



our sublayer, by contrast, viscosity is wholly negligible provided that 

A » SlJSl , (6.8) 

that is, provided that the Darcy friction from the field lines dominates the fluid friction 
from viscosity. This condition is easily satisfied, in virtue of (|6.2p . 
In this Darcy regime, p.l4p and p.lSp simplify to 

Ur = — Ku^ (6-9) 

du^ dfi dur . , 

aiid ^- -«M7r:+A^T-. (6.10) 



dz dr d 



z 



Together with p.6p and (|3.10p , equations (16. 9p and (|6.10l) describe the sublayer dynamics 
to an order of accuracy that includes the first corrections to perfect flatness. 

We now use scaling arguments, paralleling those in fj5l to verify that the corrections 
can indeed be taken as small and the sublayer treated as flat. As before, we expect that 
each term in (|6.10p is proportional to r and that 

^ = a^ + 6^r^ (6-11) 

to leading order close to the pole, where r is again dimensionless and where a^ and 6^ 
are dimensionless functions of z alone. The matching to the helium settling layer beneath 
implies that da^/dz ~ 1 and that a^ ^ ^xl^ ^^ the sublayer, if we take the constant 
value of /i above the sublayer to be zero. The condition for validity of flat, horizontally 
self-similar sublayer solutions is r^ "^ r^, say, where r^ denotes a typical magnitude of 

Within the sublayer, the pattern of dimensionless scalings (|5.5p - (|5.10p is replaced by 

Uz - S^/S (6.12) 

Ur ^ r (6.13) 

U0-rA"i (6.14) 

B'^^A'/^SJSr (6.15) 

Blr^rA^/^5^/5 (6.16) 

B'^^rA-'/'SJS (6.17) 
for all r <^ r^, r^. These dimensionless order-of-magnitude relations follow from the 
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matching to the confinement layer, again noting its horizontaUy self-similar structure, 
and from V • B' = V • u = together with (HH) and ^^, with d/dz ~ S/S^ > 1. 
Again, further detail is given in Appendix 1X1 

The horizontal velocity components in (j6.13p - (|6.14p inherit their magnitudes directly 
from those in the overlying confinement layer. For this reason, the vertical shear in the 
sublayer is larger than that in the confinement layer by a factor 6/ 6^. This shows how, 
in the limit (5^^, — > 0, preserving (J6.3I) . the sublayer regime goes over into the slip regime 
analysed in WM07. The slip regime has infinite shear at the top of the helium settling 
layer, with finite discontinuities in both horizontal velocity components. 

We can now estimate r^ and hence the flatness of the sublayer, from (j6.10p . on the 
same basis as before, namely consistency with p.lO|) and uniform downwclhng. From 
(|6.1ip we have dfj,/dr ~ bfj,r for all r ^ r^ , so that in ()6.10p the a^^ term has magnitude 
~ oifjh^r. Using ()6.12p ~ ()6.17p we find that the other terms in ()6.10p have magnitudes 
~ r((S/(5^)A^^ (the term in u^) and ^ r{5 / 5^)K (the term in Uj.). The typical magnitude 
of a^b^ can therefore be taken to be ((5/(5^) max(A, A~^). Correspondingly, since a^ ~ 
S^/S, we may define r^, the typical magnitude of a^/6^, to be af^{6^/5)'^ min(A, A~^) = 
a^(x/''7)min(A, A~^). Comparing this with (|5.1ip and recalling that a^ — N'^S^/2fl[r] 
we see that 2 

Since nx/ff ~ 0.75 x 10'^, we see that rj > r^ for a large range of N'f^ values, from 

today's value ~ N"^ down to almost three decimal orders of magnitude less. This means 

that, in most cases, our flatness assumptions hold even more strongly for the sublayer 

than for the confinement layer, in the sense of compatibility with uniform downwelling 

u^ ^u^{z) in (jXTOt . 

It is worth going beyond scale analysis to say more about the vertical structure of the 

sublayer, especially in its lower extremity or "subtail", wherein we expect |u| to decay 

exponentially with depth. Within this subtail, the helium settling layer suffers only small 

perturbations to its otherwise uniform compositional stratification d^/dz = —1. We 

denote the perturbation to /i by /i'. In the steady state, equation p.lOp may then be 

approximated as r,? / 

^ ^ /^ la^a\ 

T] OZ'' 

which can be combined with (|3.6p . (|6.9p . and ()6.10p with /i replaced by /i', to yield a 
single equation for fi' , 

a,{v/x)VJ,fi' = (A + A-i)^, (6.20) 

where V|^ = r^^d{rd/dr)/dr = V^ — d^/dz^. It is now clear that the leading-order scale 
analysis given above applies only to the main part of the sublayer and not to the subtail. 
Equation (|6.20p tells us that the vertical scale, 6i say, for the subtail must depend on the 
horizontal scale in a manner reminiscent of the heuristic boundary-layer analysis given 
in GM98. 

For instance if we assume that the actual horizontal scale is the scale r^ set by the 
confinement layer, so that a;jV|^ ^ a^r^'^, then a straightforward scale analysis of (|6.20p 
shows that, in terms of the definitions (|5.1ip and (J6.18I) . 

^)\ , (6.21) 

which is generally smaller than d^- Even at this scale, however, viscosity remains neg- 
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ligible provided that the Darcy friction from the field hnes dominates the fluid friction 
from microscopic viscosity. That is, viscosity remains negligible provided that 

A > EkV^ - Ek{d/5if, (6.22) 

equivalently Sj » (^Ik/A . (6.23) 

For realistic solar pararneters. (|6.22p is easily satisfied because Ek ^ 0.4 x 10^^. A more 
detailed analysis ( WoodI 120101 ) verifies all these properties of the subtail. The exponen- 



tially weak flow within the subtail is invisible in figure [21 and barely visible in figure [3l 
The actual horizontal scale may or may not be set by the value of r^ since, along with 
boundary conditions for (I6.20p . it will depend on the size of the downwelling region and 



indeed on how the entire global picture sketched in figure 1(b) fits together, a point to 
which we return in §111 



7. Boundary conditions for the numerical solutions 

To go beyond the self-similar, perfectly flat solutions described in [21 we must com- 
pute solutions numerically. To this end a numerical code has been written to solve the 
axisymmetric version of (|3.5[) - p.l0p in a cylindrical domain of radius r = rd, say. The 
scheme on which the code is based is summarized in Appendix JBl For reasons explained 
there and in Appendix |Cl it proves necessary to use the full time-dependent equations 
p.5|) - p.l0|) . rather than assuming magnetostrophic balance. 

To solve the equations numerically in the cylindrical domain, we need to specify bound- 
ary conditions. This inevitably involves artificial choices. The only way to avoid making 
such choices would be to fit the polar caps into the complete, and highly complicated, 
global picture. That remains a challenge for the future, requiring the quantification of tur- 
bulent processes in the tachocline and convection zone. Such quantification would have 
to include realistic descriptions of turbulent magnetic fiux pumping and of turbulent 
magnetic diffusion in the bulk of the tachocline. Also crucial is the turbulent gyroscopic 
pumping of polar downwelling and of the complementary upwelling in lower latitudes, to- 
gether with the effects on the global-scale pattern of heat flow and the resulting feedback 
on N^ distributions. 

As already noted, in the present work we are imposing a dipolar magnetic field structure 
underneath the confinement layer, and a uniform downwelling of magnitude U from a 
field-free region aloft. Field-free refers to time-averaged fields. In the example shown in 
figure m the numerical domain was defined by ^ r ^ 5, i.e. r^ = 5, and — 1 ^ z ^ 6, 
one dimensionless unit taller than shown in the figure. We imposed u^ = — 1 at z = 6 
and Br/r = 1 at z = — 1. 

As shown in gl the bulk of the confinement layer is relatively insensitive to conditions 
within the field-free region aloft, and in particular to the vertical shear du^/dz aloft. 
There, the vertical shear is tied to the temperature distribution via equation p.l6p . and 
hence to the global-scale heat flow. To avoid having to solve the complete global-scale 
problem we simply imposed "d = const, and du^/dz = at z = 6, which is consistent 
with the imposed uniform downwelling and also ensures that no viscous torque is exerted 
on the top of the domain. 

At the periphery of the domain, the artificial cylindrical surface r = Td = 5, the 
numerical algorithm requires us to impose three vertical profiles, including the vertical 
profile of Maxwell stress. The stress profile represents the field lines' connection to lower 
latitudes and the Alfvenic torque exerted therefrom. We also need to impose thermal 
and compositional stratification profiles ^(z) and ii{z) at the periphery, in a manner 
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consistent with scalings in the eonfinement layer and hehum sublayer (§^[6]). In this way 
we artificially fix the altitude of the helium sublayer. We thereby infiuence the velocity 
field as well, since it is tightly linked to the two stratifications by equations (|3.9p and 
p.lOp . The remaining peripheral boundary conditions are used to promote smoothness 
of the fields, and hence to minimise spurious boundary effects in the steady state. In 
particular, we impose ^ 



0, (7.2) 



dr\rJ 

and ^ =0 (7.3) 

or 

at r = rd = 5. 

At the bottom of the domain, the horizontal surface z = — 1, we impose the conditions 
(|3.1|) and (|3.2|) . equivalently dd/dz = 1 and d^/dz = —1. We also impose u^ = 0, i.e. 
that the interior is in solid rotation, with dimensional angular velocity $7;, as required 
by the global picture. As discussed at the end of [21 the global picture also requires 
that -B^ = 0, i.e. that there is no Alfvenic torque, on the bottom of the domain. This 
last condition cannot be directly imposed, however. Instead, it must be approached via 
iterative adjustments to the Maxwell stress profile at the periphery r = r^. The iteration 
procedure and its rationale are described in iJH 

We can now see more clearly why Uz [z) could be specified arbitrarily when constructing 
the semi-analytical solutions in 21 ^^ already mentioned, for the numerical solution we 
need to specify three vertical profiles at the periphery r = rd, which are taken to be '&{z), 
n[z), and B^{z). For the semi-analytical solutions of fJH "^i^) and ^(z) could not be 
specified independently. Rather, they were both determined by (14. 1[) . up to boundary 
conditions, as soon as Uz{z) was specified. We could still have specified B^{z), but gave 
up that freedom in order to ensure the vanishing of the Alfvenic torque as z -^ — cx), 
thus determining B^{z) via the expression (|4.9|) . Also allowed by the semi-analytical 
framework was the freedom to specify 1*0(2;) at r = rd, which we similarly gave up in 
order to ensure solid rotation as z -^ —00, thus determining u^{z) via (|4.1ip . 

More generally, within the semi-analytical framework, the peripheral and bottom pro- 
files of B^ contain equivalent information, and similarly for m^. This is because of the 
Alfvenic coupling along the field lines linking the periphery to the bottom of the do- 
main, expressed by the B • V operator in equations (|4.2|) and (|4.4|) . There is no such 
precise equivalence within the numerical framework. The time-dependence, in the equa- 
tions solved numerically, replaces static Alfvenic coupling by Alfvenic wave propagation, 
requiring one peripheral and one bottom profile to be specified, which we take to be 
Bij){z) and u^(r) respectively. This is analogous to the need for boundary conditions at 
both ends of a stretched string in motion. 



8. The numerical solutions 

Computing limitations preclude a perfect match to the real Sun's parameter values. 
They also require a slight modification to (|3.5p - (|3.10l) . explained in Appendix [B] in 
which artificial horizontal diffusivities vh, Xh are introduced. These maintain numerical 
stability while allowing small enough ly and x in the important vertical diffusion terms. 

From the scale analyses in §fj5l [6l we may identify the conditions most essential to 
reaching a qualitatively similar parameter regime — that is, qualitatively similar to a 
regime with a perfect parameter match to the real Sun. Those essential conditions are: 
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Dimensionless 


Nominal 


Value for 


Condition 


parameter 


solar value 


numerical solution 




Ro 


5 X 10"* 


10-2 


<1 


k/t] 


3 X 10* 


10^ 


>1 


a^iv/ii) 


1.4 X lO'^ 


50 


> max(A, A"^)i 


[rjr^f 


3x lO'' 


2 


> 1 


xh 


2 X 10"^ 


2 X 10-2 


< 1 


vix 


3 


1 


< 1/Ro 



\ This condition corresponds to r^ > 1, 
with r| — a^{rj/Hi) min(A, A-^) as defined in HS.lip . 

Table 1. Parameter values and conditions; see text. 



(a) The Rossby number Ro should be small in comparison with unity, so that the 
steady state is close to magnetostrophic. 

(&) The thermal diffusivity k should be large in comparison with the magnetic diffu- 
sivity ?7, so that the confinement-layer flow only weakly perturbs the background thermal 
stratification. 

(c) The confinement layer and helium sublayer should both be reasonably flat, to the 
extent that r^ and r^ are both distinctly greater than 1. With ()6.18|) in mind, we also 
take r^ > r^. 

(d) The helium diffusivity x should be small in comparison with the magnetic diffu- 
sivity 77, so that the helium sublayer is thinner than, and therefore distinct from, the 
magnetic confinement layer. 

(e) The viscosity v should be small enough that an Ekman layer does not form at 
the top of the helium settling layer, so that the flow is everywhere inviscid, even in the 
helium sublayer. With small Ro this condition is easily satisfied in the numerical model 
as well as in the real Sun, because of the factor Ro"'^ in (|6.2p . 

Leaving vh and xh aside for the moment (see Appendix |B| we can characterize the 
system by seven dimensionless parameters, including the Elsasser number A, which enters 
through the boundary conditions. Table [T] presents the other six dimensionless param- 
eters, with nominal solar values alongside the values used for the numerical solution 
presented here (figures [21 O and IH) . The last column echoes aspects of the qualitative 
parameter conditions just stated. The nominal solar values assume U = IQ-^cms"^. 

In order to allow the stratification surfaces to develop a significant tilt, the horizontal 
size of the numerical domain, rj, was chosen to be of the same order as r^)] specifically, 
we chose r^ — 5. However, the precise value of r^, as determined by (|5.1ip . depends on 
the precise value of A, which is set indirectly via the boundary condition for Br- For 
the case shown in figures [2l El and El this boundary condition was Br/r = 1 at z — —1, 
and we find from the solution that A = B'f\ k, 3.5. With the parameter values in the 
second- last column of Table [1] this in turn means that r^ w 50/3.5 « 14, so in fact r^} is 
slightly smaller than r^. Nevertheless the numerical solution is remarkably close to being 
flat, even close to the periphery of the domain r = r^. Some effects of departures from 
flatness can be seen near the periphery in, for instance, figure (TUl below, but they do not 
qualitatively alter the nature of the flow. The conflnement-layer dynamics can therefore 
apply, at least qualitatively, even at colatitudes for which the semi-analytical solutions 
of ijH are not strictly valid. We should therefore regard (j5.12|) as a conservative estimate 
of the range of interior field strengths for which the confinement-layer regime applies. 

The vertical profiles of d, jjl and B^ at the periphery of the computational domain were 
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Figure 5. A numerical solution of the confinement-layer equations with r| ~ 14 and Elsasser 
number A « 3.5. Other parameter values are given in the second-last column of Table [l] From 
the same solution as shown in figure [2] Both here and in figure [2l a top layer 5 ^ z !^ 6 has been 
omitted from the plots. Figures (2] [5] and [6] are all views of the same numerical solution. 

initially taken from a semi-analytical solution. The resulting steady-state meridional flow, 
and the distribution of Coriolis torque it exerts on each field line, cannot be precisely 
known in advance. So the steady state found, with this choice of the peripheral B^ profile, 
will generally be such that some of the balancing Alfvenic torque is exerted on the bottom 
of the computational domain, rather than on the periphery r — r^i. 

As mentioned in the previous section, we can eliminate this bottom torque through 
iterative adjustments to the torque at the periphery. Experience with the numerical solu- 
tions reveals that such adjustments are propagated along the magnetic field lines without, 
in most cases of interest, greatly perturbing the poloidal velocity and magnetic field com- 
ponents. So we can reduce the bottom torque simply by mapping or transferring that 
torque from the bottom to the periphery of the domain, along the field lines. The bottom 
torque is thereby reduced, though not all the way to zero since the other fields adjust, 
slightly reshaping the field lines. The process can then be repeated, further reducing the 
bottom torque. After a sufficient number of iterations, the bottom torque can be reduced 
to negligible values, with practically all the torque transferred to the periphery. 

Figure [5] shows plots of the steady-state streamlines and magnetic field lines from the 
numerical solution whose parameter values are listed in Table [l] and whose meridional 
cross-section was presented in figure [2j The bottom torque is negligibly small. Figure |6] 
shows the vertical profiles of Uz, u^/r, Bz and B^/r on the rotation axis, from the same 
numerical solution. 



9. Upper-domain "slipperiness" 

On the rotation axis, where the profiles in figure|6]were taken, the stratification surfaces 
are flat for any finite a^ and a^j . If the numerical solution were in perfect magnetostroph- 
ic balance then we could use its Uz{z) profile to calculate the other field components on 
the axis by the procedure for constructing semi-analytical solutions described in §4] But 
the numerical solutions are not in perfect balance, especially toward the upper part of 
the domain, where the Lorentz and Coriolis forces become small and the artificial viscous 
forces relatively more significant, along with numerical truncation errors and other small 
effects. In particular, the numerical Uz{z) and Ur{z) profiles will not conform precisely 
to the decay law (I4.10|) as z increases. So the semi-analytical solution obtained by this 
process cannot perfectly match the numerical solution, even on the rotation axis. Indeed, 
such a semi- analytical solution will often exhibit wild deviations from the numerical 
solution toward the upper part of the domain. There, the delicate balance of terms 
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Figure 6. Vertical profiles from a numerical solution of the confinement-layer equations with 
r^ ~ 14 and Elsasser number A ~ 3.5. Other parameter values are given in the second-leist 
column of Table [l] From the same solution as shown in figure (2] 

gives the dynamics a certain "slipperiness" , as already evidenced by the upper-domain 
sensitivity of u^ and du^jdz to values of the decay constant 7 in ()4.10|) . 

To enable a meaningful comparison between the numerical and semi-analytical so- 
lutions, we are therefore compelled to make small adjustments to 1*2(2) in the upper 
domain, to make u^iz) and u^i^z) conform to (|4.10p . before using them to compute a 
self-similar solution. In the case shown here the decay constant 7 was chosen, purely for 
illustration, to be 2.24. The required adjustment to Uziz) is then very small indeed; the 
solid Uz curves on the left of figures [3] and [6] are practically indistinguishable from each 
other. 

As already mentioned, conditions aloft in the real Sun may well involve small-scale 
MHD turbulence. The resulting departures from magnetostrophic balance may well pro- 
duce asymptotic behaviour aloft that disagrees with all the solutions obtained here, 
whether semi-analytical or numerical. A realistic matching to conditions aloft remains a 
challenge for future modelling work. 



10. Confinement layers with no helium settling layer 

Although the presence of the helium settling layer influences the structure of the con- 
finement layer in today's Sun, it is not actually essential to magnetic confinement. The 
qualitative picture sketched in figure 1(b) might therefore apply also to the early Sun, 
before the helium settling layer developed. 

The analysis presented in ij4]holds good for any suitable profile of u^iz), provided only 
that the thermal stratification is sufficiently strong. Compositional stratification enters 
only indirectly, via the shape of u^i^z) in its lower part representing the helium sublayer 
and subtail. As pointed out below (|6.2ip - (|6.23|) . the shape depends in turn on how the 



confinement layer and sublayer fit into the global picture sketched in figure 1(b) 



For the early Sun, with no compositional stratification, we have a similar indeterminacy 
in the shape oiuz{z)- However, we may note that with no helium settling layer the scaling 
analysis in ij5]then applies not only within the confinement layer, but also in the region 
immediately beneath. An argument similar to that leading to the estimate bi of the 
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Figure 7. An semi-analytical solution of the confinement-layer equations with no helium sett- 
ling layer. The downwelling was chosen such that j = 2.24. The u^ profile has been rescaled by 
a factor of 10. The z-origin is arbitrary. 




Figure 8. The magnetic confinement layer near the north pole in a model for the early Sun, 
with no helium settling layer. The plot is from the same semi-analytical solution as that of 
figure [T] 



subtail scale in (|6.2ip predicts that |u| decays with depth on a lengthscale ^ S, rather 
than Sf. 

Analytical solutions with no helium settling layer can be calculated in the same way 
as in [21 Figure [7] presents such a solution. In this case, the downwelling profile was 
chosen to be Uz{z) = —(1 -|- exp(— 7(2 — 2)))^-^/''', so that Uz = 0(exp(2;)) as 2; ^ —00 
and Uz = —1 + 0(exp(— 7z)) as z -^ +00. We have taken 7 — 2.24, again purely for 
illustration, to allow a more direct comparison with figures [3] and [51 Figure [5] shows a 
vertical cross-section through the solution presented in figure [71 
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11. Conclusions and future directions 

We cannot yet claim to have a complete tachocline theory. Indeed, the confinement 
layer and helium sublayer form only two pieces of a complicated jigsaw puzzle. Other 
aspects of that jigsaw include the way in which the confinement layer matches upward 
to the relatively large negative shear in the bulk of the tachocline, and the way in which 
the baroclinic temperature anomalies induced by the tachocline's MMC fit into the per- 
turbed global-scale heat flow. In particular, without putting the whole jigsaw together 
we cannot quantitatively predict the thickness of the tachocline. Nor can we predict the 
precise shapes of the vertical profiles of B and u in the confinement layer. Those pro- 
file shapes depend on matching to conditions not only aloft but also equatorward, where 
stratification surfaces and field lines extend into colatitudes outside the polar downwelling 
regions. 

However, the results obtained here give us the first fully consistent model of polar field 
confinement, as such, together with insight into how it could work in today's Sun. The 
results cover a wide range of possible downwelling values and interior field strengths (end 
of ©. We have also shown, in flOl how confinement could have worked in the early Sun. 
The dynamics is similar apart from the slightly deeper penetration of the MMC in the 
absence of the helium settling layer and helium sublayer. We can use the resulting in- 
sights, alongside our well-established understanding of the gyroscopic pumping of MMCs, 
to say something new about the early Sun and the solar lithium-burning problem. 

Standard solar-evolution mo dels predict surface li thium abundances higher than ob- 
served by a factor ^ 10^ (e.g. IVauclair et aLlll978[) . The reason is that the standard 



models mix material down to the bottom of the convection zone but no further. To 
destroy lithium, material from the convection zone must be mixed or circulated to some- 
what greater depths and therefore to somewhat higher temperatures, beyond those at the 
bottom of today's tachocline. However, there is no evidence of depletion of the convec- 
tion zone's beryllium, which is destroyed at only moderat ely higher temperatures than 



lithium . Furth er discussion and references may be found in IChristensen-Dalsgaard et 



(I994), and in IWoodI (2010!, chap. 6). Here we argue that a quantitative version of the 



scenario sketched in figure 1(b) has promise as a way of circulating material to the re- 
quired depth in the early Sun, and no further, thus making sense of the high beryllium 
as well as the low lithium abundance. 

As already mentioned in fjl] the downwelling MMC in the polar tachocline that makes 
field confinement possible can be regarded as due to a gyroscopically-pumped MMC 
trying to burrow downward, but held in check by its encounter with the interior field Bi 
and with the helium settling layer, if present. If, in a thought experiment, we were to 
switch off the interior field Bj, then the downwelling would spread or burrow to ever- 
increasing dept hs. The timcscalc for such burrowing is inversely proportional to il? (e.g. 



Mclntvrdl2007l equation (8.15)ff.); one may think of rotational stiffness as strengthening 
the burrowing tendency. 

Now, because the early Sun rotated much faster than today, not only would there 
have been no helium settling layer but, also, the burrowing tendency would have been 
much stronger than today, tending to push the bottom of the tachocline downward. This 
reopens the possibility conjectured in GM98 that there might have been a ventilated 
"polar pit" in which most of the convection zone's lithium, though not too much of 
its beryllium, was burnt during the first gigayear or two of the Sun's main-sequence 
evolution. 

To take this further we again need to consider the way in which the confinement layer 
fits into the global picture. It is arguable that the bottom of the entire polar downwelling 
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Figure 9. Schematic drawing of the magnetic confinement layer and its immediate surroundings 
at the bottom of the high-latitude tachocline. Close to the pole the interior magnetic field (solid 
lines) is confined by the downwelling MMC (dashed streamlines). The vertical scale has been 
greatly exaggerated. 



region is depressed relative to its surroundings, forming not so much a "pit" as a shallow 
"frying pan" , too shallovif to burn lithium in today's Sun but possibly just deep enough 
in the early Sun. 

Here we need to distinguish the shape of the ventilated region from the shapes of 
the stratification surfaces, which latter must remain relatively flat, meaning close to the 
horizontal. Figure [9] sketches the way in which the confinement layer might fit into its sur- 
roundings near the bottom of the polar tachocline. The stratification surfaces are shown 
dotted. At the periphery of the polar downwelling region, the field lines (solid) spiral 
outward and upward from the confinement layer on their way to lower latitudes. They 
will tend to splay out, as well as slanting upward, as they emerge from the downwelling 
region. The MMC will similarly slant upward, flowing approximately along the field lines 
(dashed streamlines). This is because the splaying-out increases the magnetic Reynolds 
number beyond the order-unity values characteristic of the confinement layer. Further 
out, the field lines must continue to rise through the tachocline until they encounter the 
convection zone's overshoot layer, where they are held horizontal by turbulent magnetic 



flux pumping as suggested in figure 1(b) On the way we must expect turbulent eddy 
fluxes to become increasingly important, decoupling the MMC's upwelling streamlines 
from the time-averaged field lines and leaving the upwelling free to spread over a wide 
range of latitudes, constrained only by mass conservation and global-scale heat flow. 

Such a picture applies equally well to today's Sun and to the early Sun, the main 
difference being that the ventilated polar region (unshaded in figure E]) is likely to have 
been pushed deeper in the early Sun with its much faster rotation, stronger burrowing 
tendency, and global-scale |Bi| values only modestly larger. The ventilated polar regions 
could well have been deeper by many tens of megametres, as required to burn lithium|jj 
The early Sun would have started to form its helium settling layer just below these 
ventilated, lithium-destroying polar regions, i.e. just below the polar confinement layers. 
Then, with the gradual diminution of ilj through solar-wind braking, the confinement 
layers, marking the bottom of the ventilated regions, would have retreated upward, and 
the top of the helium settling layer would have followed them upward as new helium 
strata formed. 

Within the peripheral lightly-shaded region in figure |9l into which the MMC does 
not penetrate, we suggest that ventilation is weak or nonexistent and that shear will be 



f This deepening is additional to the deepening of the convection zone itself, in the early Sun 
relative to today' s Sun, amounting to several more tens of megametres according to standard 
solar models (e.g. ICiacio et a/.lll997f ). 



Confinement of the Sun's interior magnetic field 29 

limited by the Ferraro constraint. The darker shading represents the top part of today's 
hchum settling layer. 

As the lightly-shaded region expands upward and outward beyond the immediate sur- 
roundings sketched in figurelHl through the tachocline towar d the overshoot layer, we may 



surm i se that small-scale MHD instabilities will kick in fe.g. lSpruitJl2002tlGilman fc Callv 



2007t IParfrev fc Menoull2007l . & refs.), breaking the Ferraro constraint and blurring the 



distinction between the shaded and unshaded regions as turbulent eddy fluxes increase. 
So a larger-scale picture of the "lithium frying pan" would show its upward-sloping lower 
boundary becoming increasingly porous and indistinct at greater colatitudes. 

The global tachocline model that would be needed to test, and to begin to quantify, 
the foregoing speculations would have to describe 

(a) the precise way in which turbulent stresses in the convection zone and tachocline 
gyroscopically pump the polar downwelling needed to confine B; in polar latitudes; 

(6) the global-scale distribution of temperature and heat flow that fits in with the 
MMCs; 

(c) the turbulent magnetic flux pumping by convective overshoot assumed to confine 
Bi in extra-polar latitudes; 

(d) the extent to which the winding-up of the time-averaged toroidal field in extra- 



polar latitudes (figure 1(b)) is hmited by turbulent eddy fluxes; 

(e) the reaction of the overlying turbulent layers to all of the above, especially the 
deficit in the convection zone's differential rotation governing the torques exerted from 
above, for instance via feedback on the strength of gyroscopic pumping of the MMC. 

Progress on these formidable problems will depend on finding suitable ways to model 
the turbulent processes. 
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Appendix A. HELIUM SUBLAYER SCALINGS 

For perturbations B' to a background magnetic field Bj, the Lorentz force may be 
written as 

F = (Vx(Bi + B'))x(Bi+B') = -V(i|B,|2 + BrB' + i|Bf )-H(Bi-fB')-V(Bi + B') 

(Al) 
If, as here, the background field Bj is curl-free, then all the terms quadratic in B; vanish 
and we have 

F = -V(Bi • B' + i|B'p) + B' ■ VBi + (B; + B') • VB' . (A 2) 

If Bi and B' are axisymmetric, and have the scalings given by (I6.15|) - (|6.17p . with d/dz ^ 
S/6^ ^ 1, then all components of B'- VB' are smaller than the corresponding components 
of Bj • VB' by factors {5^/5)'^, with one exception. The r component of B' • VB' includes 
a term B'^/r, the divergence of the hoop stress. Relative to the r component of Bi • VB' 
this term is of order {5^/5)'^ /h? . However, the hoop-stress term itself is smaller than the 
r component of the Coriolis force by a factor {S^/6)'^, even for small A, because of ()6.14|) . 
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We also find that, because Bi^ — 0, the poloidal components of V(Bi ■ B') exceed in 
magnitude those of V(5|B'p) by a factor (6/6^)^. Their azimuthal components are both 
zero since we consider only axisymmetric fields here. We may therefore neglect all terms 
quadratic in B', in the asymptotic regime (j6.3p . so that (|3.13p becomes 

e^ X u = -Vp' + a^de,. ^ a^/ie^ - V(Bi • B') + B' • VB; + Bi • VB' . (A 3) 

The flow through the sublayer produces only a small perturbation to the otherwise uni- 
form thermal stratification. From p.9|) and p.lOp we see that, within the sublayer, varia- 
tions in d are smaller than variations in ^ by a factor x/k ^ 1. The thermal stratification 
in the sublayer may therefore be treated as horizontally uniform, allowing the thermal 
buoyancy term a^ "de^ in (jA3|) to be incorporated into the pressure field along with the 
— V(Bi ■ B') term. We denote the modified pressure field as p. 
The perturbed steady-state induction equation is 

= (Bi -f B') • Vu - u • V(Bi + B') + V^B' , (A 4) 

again on the assumption that Bi is curl-free. It is readily verified from the scalings ()6.12p - 
(|6.17p and d/dz ~ 5/5^ ^ 1 that each component of u • VB' is of the same order as 
the corresponding component of B' • Vu, but smaller than Bi • Vu by a factor [5^/5)'^ . 
Furthermore, provided that the horizontal scales r^ and r^ are both ^ 1 we may also 
make the boundary-layer approximation, V^ ~ d"^ jdz^ . Thus, in the asymptotic regime 
(|6.3|) . we may simplify (|A4[) to 

= Bi- Vu-u- VBi + — ^B' . (A 5) 

If the field Bi were uniform, and directed along the axis of rotation, then (jA3p and 
(|A5|) would reduce immediately to (|6.5p and ()6.4p respectively. Since Bi is axisymmetric 
and smooth, this reduction still holds as a first approximation within some neighbourhood 
of the axis. It is sufficient to show that this neighbourhood includes the entire sublayer 
within a radius r ~ r^. We first note that, since V • Bi = 0, we have B{r/r = —^dBiz/dz 
as in (|4.8I) . In the sublayer we have dBiz/dz ^ Biz (dimensionally, dBiz/dz ^ Biz/5) as a 
consequence of the matching to the confinement layer, as can be seen, for instance, from 
the solid curve in the right-hand panel of figure[31 Now applying the scalings (|6.12p - (|6.17p 
we find that, in (jA3p . each component of B' ■ VBi is smaller than the corresponding 
component of Bi • VB' by a factor S^/5, and that, in (jASp . each component of u • VBi is 
smaller than the corresponding component of Bi • Vu by the same factor 6-^/d. Moreover, 
even at colatitudes r r^ r^j, the contributions Bizd'B' /dz and Bizdu/dz dominate all other 
contributions to Bi ■ VB' and Bi • Vu by factors of at least 6/6^. So (|A3p and (lASp do 
indeed reduce to (|6.5p and (|6.4I) . 

For r <r^, tilting of the compositional isopleths produces variations in the dimension- 
less altitude of the sublayer no greater than 0{6^/6), so that Biz may be assumed constant 
within the sublayer as assumed in the derivation of (|6.9I) and (j6.10p . Finally, we note that 
the foregoing picture applies not only to steady states but also to time-dependent states 
with any timescale, such as 6/U = S'^/rj, that is long in comparison with the timescale 
S'i/ri for magnetic diffusion across the sublayer. On any such timescale the sublayer 
therefore behaves like a porous medium. 



Appendix B. THE NUMERICAL SCHEME 

We wish to solve a suitable version of equations p.Sp - p.lOp in axisymmetric cylindrical 
polar coordinates. We introduce streamfunctions ^f and A, i.e. azimuthal vector-potential 
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(Bl) 

(B2) 



guaranteeing that the fields are divergence- free; r'^ is sometimes called the Stokes stream- 
function. The azimuthal vorticity w^ and electric current density J^ = (V x B)^ are 
related to ^I' and A by 



and 



W0 = 



r-2)^ 



N^-r-^)A 
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As already explained, and further discussed below, we introduce anisotropic viscosity 
and helium diffusivity with dimensionless horizontal components vh and xh- So equa- 
tions p.5p - p.lO|) are replaced by 



^^iy{ru,) 



Ro 



Dt dz 

T){io^/r 



(9* 1 

— = -B-V{rB^) 



d{ru^,u^/r) 
d{z,r) 



-a,', 



f Ek 

_ du^ 

dz 

d^i 

' dr 



d^ 

dz^ 



VH 






U0 



(B5) 



rB • V{J^/r) 



d{rB^,B^/r) 
d{z,r) 



Ek 



dz 



VH V 



m 

1 D(rA) 
r T>t 



= rB ■ V{u^lr) 



A 



H 



V'- 



W0 



Bd 



92 



D/i 



dz 



XH^Ji 



M 



(B6) 
(B7) 
(B8) 
(B9) 
(BIO) 



We solve these equations using a simple finite-difference scheme on an Eulerian grid regu- 
larly spaced in r and z at intervals Ar and Az. The outer boundary of the computational 
domain is at r = rj. The inner boundary is at r = 2Ar, i.e. two grid intervals from the 
coordinate singularity at the rotation axis. Because of the directionality of operators like 
u • V and B • V, the spatial derivatives are calculated using two-point, one-sided (first- 
order) finite differences whose directions are chosen to ensure numerical stability at the 
grid scale. For reasons of symmetry and good behaviour near the coordinate singularity, 
the finite differencing is done by locally approximating the fields ^/r, u^/r, uj,f,/r, A/r, 
B^/r, 1}, and /i as functions that are linear in z and in r^ over a single grid interval. This 
ensures that the error is 0(Ar) even for small r. Field values for r < 2Ar are obtained by 
extrapolation from r = 3Ar and r = 2Ar, again assuming linear functional dependence 
on r^. 

With the parameter values given in Table [1] the dimensionless helium-sublayer and 
Ekman-layer thicknesses are S^/S = {xIvY^"^ ~ 0-14 and 6^]^/5 — Ek^" w 0.01 respec- 
tively. We have chosen a vertical grid interval Az — 0.01, dimensionally 0.01(5, which is 
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small enough to resolve the helium sublayer accurately. This Az is too large to resolve 
any Ekman layers. However, Ekman layers are prevented from becoming significant by 
careful choice of the code representing the boundary conditions. By allowing slip veloci- 
ties and making viscous stresses negligible at the boundaries, we have been able to keep 
Ekman-layer formation so weak as to play no significant role in the dynamics. Uniform 
rotation is imposed at and just beneath the bottom of the domain, via the B • V term 
in the azimuthal component of the induction equation (|B 7p . 

An explicit Eulerian timestepping scheme is used to evolve the system. The timestep 
At must be small enough to resolve thermal diffusion at the grid scale (which is the fastest 
process at this scale and therefore determines the Courant-Friedrichs-Lewy condition). 
So from TablelH At < {ri/Hi){Azf = IQ-^ x (O.Ol)^ = 10-^ dimensionally IQ-^S/U or 
10-4(2^)"^. 

We use a semi-analytical confinement-layer solution, of the kind described in iJH as 
the initial condition. The system typically takes several domain-scale magnetic diffusion 
times to reach a steady state. As explained in fJSl multiple iterations of the peripheral 
Brf,{z) profile are then required to achieve a steady state with vanishing B^ir) at the 
bottom. To make the computation feasible, in a domain wide enough to accommodate 
noticeable tilting of the stratification surfaces, we have used r^ — 5, dimensionally 5(5, 
and a horizontal grid interval Ar = 0.1, dimensionally 0.1(5, larger than the vertical 
grid interval Az by a factor of 10. For numerical stability, the dimensionless horizontal 
viscosity vjj and helium diffusivity xh rnust then be chosen so that the diffusive terms 
in (jB 5[) . (IB 61) and (|B 10[) dominate the advective terms at the gridscale. In practice, 
we found it sufficient to take vh — Xh — 10, i.e. to make them a factor 10 greater 
than the corresponding vertical diffusivities, as suggested by advective-diffusive scaling 
when Ar / Az — 10. We have verified, in smaller computational domains, that the coarser 
horizontal resolution and anisotropic diffusivities do not qualitatively affect the steady 
state of the system. 

At each timestep, the azimuthal vorticity oj^ is updated and the streamfunction ^ 
then computed from (|B 3p by inverting the operator V^ — r~^, approximated using cen- 
tred differences. The inversion is perfor med iteratively, using a successive-overrelaxation 



method described in lPress et al\ ( 19861 ). During the early evolution, when the dynamics 
is dominated by timescales not much longer than the timestep At, many such iterations 
are required, at each timestep, to achieve convergence. At later times the same degree of 
convergence can be achieved with far fewer iterations. Since we are interested only in the 
ultimate steady state, we can tolerate a larger error in the inversion durin g the system 's 



transient evolution. Further details of the numerical code are spelt out in IWoodI ( 2010r ). 



As anticipated from the small values of Ro and Ek, the steady state is found to be 
close to magnetostrophic balance: the motion is scarcely distinguishable from one in 
which the balance conditions (|3.14p and (|3.15l) hold exactly. Indeed, in p.lSp the terms 
in a^ and a^ dominate so strongly that, in the case of figure [2] for instance, the tilting is 
barely visible, as has been verified from plots, not shown, of the thermal as well as the 
compositional stratification surfaces. This of course is no more than was anticipated from 
the scaling arguments of §33111 So the balance p.lSp holds in an almost trivial sense, 
with the remaining terms scarcely able to produce any noticeable effect. The azimuthal 
momentum balance described by p.l4p is less trivial, but here also we find that the 
left-hand side and right-hand side are close to being equal, in the present notation 

19* 1 

--^«-B.V(ri3^). (Bll) 

r oz r^ 

Figure [10] shows the left-hand side and right-hand side of (|B lip for the case of figure [21 
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Figure 10. The left and right panels show contours of the left-hand side and right-hand side 
of (|B 11|) respectively, in a vertical cross-section for the numerical solution of figures [2l [5] and 
El The bulk of the confinement layer is close to magnetostrophic balance. The largest fractional 
departures from balance occur in the regions above and below the confinement layer where both 
sides of (|B 11|) are numerically small. Unlike that in figure [21 the cross-sections shown here span 
the entire numerical domain including the top part 5 ^ « ^ 6. 



The balance is quite accurate despite the modest Ro value chosen for the numerical 
solution, 10^2 i^^ Table [Tj 

It might b e thought tha t imposing magnetostrophic balance from the start, as first 
suggested by iTaylon (|l963l ) , would filter out all the fast oscillations — including inertial 
or epicyclic oscillations as well as Alfven waves, gravity waves, and the various hybrid 
types — and thereby allow larger timesteps to be used. However, th e imposition of mag - 
netostrophic balance leads to pathological behaviour at small scales ( Walker et al\\l99&i] . 
Far from eliminating or slowing the fast oscillations, the imposition of balance exacerbates 
the problem, for reasons explained in Appendix [Cl 



Appendix C. MAGNETOSTROPHIC BALANCE AND 
NUMERICAL ILL-CONDITIONEDNESS 

Filtering out fast oscillations by imposing some kind of balance is a familiar, and 
often effective, device in many other problems involving stiff differential equations. A 
well known example is that of fiuid flow in non-MHD fiuid systems with strong rotation 
(small Ro) and stable stratification. The standard "quasi-geostrophic equations" result 
from imposing geostrophic or Coriolis-pressure-gradient as well as hydrostatic balance, 
thereby filtering out inertia and gravity waves as well as sound waves. The filtering saves 
computational resources by allowing relatively large time steps. 

Such filtering turns out, however, to be ineffective in the confinement-layer problem. 
Indeed — at first sight paradoxically — the imposition of magnetostrophic balance ex- 
acerbates the timestepping problem. Far from eliminating fast o scillations, it introdu ces 
spurious modes of oscillation that are even faster, as s hown by Walker et all (|l998[) in 
the context of the terrestrial dynamo problem. Following I Walker et all we show how this 
pathology can be understood through an idealized analysis of the fast oscillations, first 
in the unfiltered and then in the filtered equations. 

The reason for the pathology is the interplay between the Coriolis and Lorentz forces. 
Stratification 7V^ is relatively unimportant, as will be shown shortly. We therefore start 
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with the hnear theory of MC (magneto-Coriohs) waves, i.e. smaU plane- wave disturbances 
to an unstratified, incompressible fluid with solid rotation fl and a uniform magnetic 
field B. Neglecting viscosity and magnetic diffusivity, we find the well-known dispersion 
relation 

w2-(B-k)2 = ±2rj-ka;/|k|, (CI) 

where w is the frequency and k is the wavevector, both dimensional here. If we take 
the limit of rapid rotation, |ri| — > oo, then for most choices of k the four roots of this 
dispersion relation are asymptotically 

and ,,^±(|J^|k|. (C3) 



The modes corresponding to (|C 2[) are inertial or epicyclic waves — in this context some- 
times called "fast MC waves" — and those corresponding to (|C3p are "slow MC waves" . 
By imposing magnetostrophic balance we neglect relative fluid accelerations, which cor- 
responds to dropping the w^ term from the left-hand side of (IC 1|) . The dispersion relation 
then becomes /-p, , n2 

-■-^IttM^ (C4, 



So imposing magnetostrophic balance eliminates the two "fast" branches (jC 2p of the full 
dispersion relation (jC 1[) . 

However, not all modes of the full dispersion relation (|C ip have the asymptotic be- 
haviour described by (|C 2|) and (|C 3p . Even in the presence of rapid rotation, there are 
always some modes whose k values satisfy 

|B-k| » |2n-k|/|k| (C5) 

by an arbitrarily large factor. For instance we can fix the direction of k and make |k[ 
arbitrarily large. Such modes behave like Alfven waves, with uj « ±B • k. Imposing 
magnetostrophic balance removes the mechanism for Alfven wave propagation, and must 
therefore alter the behaviour of these modes. In fact their frequencies become arbitrarily 
higher than Alfven wave frequencies. This can be seen at once by inspection of (jC4l) and 
(|C 5[) . In summary, even in a rapidly rotating system some modes of the full dispersion 
relation always feel the Lorentz force more strongly than the Coriolis force, i.e., they 
satisfy (|C 51) . and these modes become ill-behaved under the assumption of magnetostroph- 
ic balance. A numerical scheme that imposes magnetostrophic balance while retaining 
realistic (Laplacian) magnetic dissipation will therefore be ill-conditioned. 
If we introduce stratification N'^ then (IC4I) becomes 



-=±^[(B•kf|kp+7v^lk^^^/^ (C6) 

where k/f is the horizontal projection of k. Therefore the stratification (a) makes little 
difference to the large-|k| behaviour but (b) always increases the frequency of the ill- 
behaved modes and thereby, if anything, further exacerbates the problem. 
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